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1 Introduction 

In this paper, we report on an implementation in the free software Mathemagix [9! of 
lacunary factorization algorithms, distributed as a library called LacunaryxQ These 
algorithms take as input a polynomial in sparse representation, that is as a list of 
nonzero monomials, and an integer d, and compute its degree-d factors^ The com¬ 
plexity of these algorithms is polynomial in the sparse size of the input polynomial 
and d: The sparse size of a polynomial / = jjj =1 • • • X* n,i E Q[Xi,..., X n ] is 

Lj{h{cj) + E«log(l + a i,j)) where the /-;(■) denotes the height of a rational number, that 
is the maximum of the absolute values of its numerator and denominator. 

In our implementation, we focus on polynomials with rational or integer coef¬ 
ficients. Algorithms in the univariate case have been given by Cucker, Koiran and 
Smale [3J and generalized by Lenstra j8). In the multivariate case, the first algorithms 
were given by Kaltofen and Koiran |6][7l. We described simpler algorithms for mul¬ 
tilinear factors, with a different approach, with Chattopadhyay, Koiran, Portier and 
Strozecki |2j[iJ and then generalized them to bounded-degree factors I4 5j. We imple¬ 
ment a variant of Lenstra's algorithm for the univariate case, and our most recent algo¬ 
rithms for the multivariate case. The algorithms of Koiran and Kaltofen are not easily 
implementable since they rely on the value on some non-explicit number-theoretic 
constant. 

*Partially supported by a LIX-Qualcomm®-Carnot postdoctoral fellowship, and by the French 
Agence Nationale de la Recherche under grant CATREL #ANR-i2-BSo2-ooi. 

I This library and its documentation can be found at http: //mathemagix. org/www/lacunaryx/doc/html/index. en. html 

2 Here and henceforth, the expression degree-d factor means an irreducible factor of degree at most d. 
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2 Algorithms 

This section is devoted to the description of the algorithms we implemented. We 
refer to the original papers for more formal descriptions, proofs of correctness and 
complexity estimates. 

2.1 Lenstra's algorithm 

Lenstra's algorithm has two disjoint steps for computing cyclotomic factors on the one 
hand, that is factors that divide X r — 1 for some r, and non-cyclotomic factors on the 
other hand. Let us first describe the step for non-cyclotomic factors. 

Theorem i (Lenstra's Gap Theorem) Let d be a positive integer arid f = f 1 +f 2 e Q[X], 
There exists a constant 7(/, d ) such that ifwdl^ff) — deg(/i) > 7(/, d), every non-cyclotomic 
factor of degree at most d of f divides both f\ and f2- 

The constant in the theorem is polynomial in the sparse size of / and in d. It yields 
an algorithm to compute the degree-d non-cyclotomic factors of a given / 6 Q[X] in 
time polynomial in the sparse size of / and d: 1. Compute y(f, d) and f\ of minimal 
degree such that val(/— ff) — deg(/i) > 7(/,d); 2. Recursively compute /2, ..., 
f s such that / = fi + ■ ■ ■ + f s with the same property; 3. Compute the factors of 
gcd(/i,.. .,/s). The theorem ensures that every degree-d non-cyclotomic factor of / 
is a common factor of /1, ..., / s , that is a factor of gcd(/i,... ,/ s ). To compute the 
multiplicities, the same algorithm is applied to the derivatives of /, or more precisely 
to its sparse derivatives: The sparse derivative /-’- of / is defined as (//X val ^) ; where 
' denotes the standard derivative. 

The computation of the cyclotomic degree-d factors of / works as follows: If g is 
cyclotomic, it divides by definition X r — 1 for some r. The cyclotomic factors of / that 
divide a fixed X r — 1 also divide gcd (/, X'' — 1 ). Let us write / = Yy =1 CjX a i and define 
/mod r = £* =lC/ X*; mod L p hen jr mod r sat i s fi e s gcd (/, X r - 1 ) = gcd(/ mod r , X r - 1 ). 
This shows that for each r, one can compute the common factors of / and X' — 1 in 
polynomial time in r and the sparse size of /. Number-theoretic considerations show 
that all degree-d cyclotomic polynomials divide some X r — 1 with r < 2 d 2 . Altogether, 
the degree-d cyclotomic factors of / can be computed in polynomial time in d and the 
sparse size of /. Again, the multiplicities are computed using the same algorithm on 
the sparse derivatives of /. 

2.2 Variant of Lenstra's algorithm 

To describe our variant of Lenstra's algorithm, let us reformulate it slightly. With 
the notation of the previous section, let g = gcd(/i,... ,/ s ), and h = f/g. Lenstra's 
algorithm can be described as follows: l. Compute g and h; 2. Factor g; 3. Compute 
the cyclotomic factors of h; 4. Return the union of these factors (viewed as multisets, to 
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take multiplicities into account); 5. Apply the same algorithm to the sparse derivatives 
of/. 

We build on this view in our algorithm. Instead of computing a single g and a 
single h, we actually compute two sets G and H such that for any degree-d polynomial 


l, 


mult ff) = 1 £geG mul ttC?) + LheH™\t e (h) 

\Eg6G mul t£(g) 


if i is cyclotomic, and 
otherwise 


(1) 


where mult//) denotes the multiplicity of £ as a factor of /. Let us remark that 
the most expensive part of Lenstra's algorithm is the factorization of the polynomial 
gcd(/i,... ,/s). Our modification allows us to reduce as much as possible the degree 
of the polynomials that have to be eventually factored. 

To compute the sets G and H, we reverse the logic of Lenstra's algorithm: While 
Lenstra computes an a priori bound 7 (/, d) to split / as a sum, we split / as much 
as possible (say as a sum of monomials) and then merge together the summands, 
beginning with the closest ones. At each step, if the summands have a non trivial 
gcd g, we apply recursively our algorithm to // g to get the two sets G and H, and 
return (G U {g},H). Else, we merge the two closest summands (that is we replace 
both of them by their sum) and compute again the gcd of the new summands, until 
there remains a unique summand. As described, the algorithm would always return 
an empty set H since the gcd of the summands is computed at each step. Actually, 
this is useless in the following case: If two summands /• and / !+ i have been merged 
but val(/j + i) — deg(/■) > 7(/■ + / + j, d), Lenstra's Gap Theorem ensures that the non- 
cyclotomic degree-d factors of /■ + f l+ \ are common to / and f 1+ j. The gcd does not 
need to be computed, and we can directly merge the next summands. When it remains 
a unique summand, it goes into H rather than G if the gap that has been merged is 
large enough. The formal description of this algorithm is given as Algorithm [1] Note 
that the sets G and H give a partial factorization of / in the sense that / = YI(cc,jh 

After G and H have been computed, it remains to compute all the factors of the 
polynomials in G using any classical algorithm, and the cyclotomic ones for polyno¬ 
mials in H using a variant of Lenstra's algorithm for cyclotomic factors. Finally, the 
union of these factors is returned. We remark that our method directly provides the 
multiplicities of the non-cyclotomic factors. 

Let us describe our variant of Lenstra's algorithm for cyclotomic factors. Let h 6 H 
be a polynomial, and suppose we aim to compute its degree-d cyclotomic factors. 
As mentioned in the previous section, a cyclotomic factor divides X r — 1 for some 
r. Lenstra proved that r < 2 d 2 for a degree-d cyclotomic factor. Actually, one can 
compute a much smaller set R C { 1 ,..., 2 d 2 } such that each degree-d cyclotomic poly¬ 
nomial divides X r — 1 for some r E R. By definition, the r-th cyclotomic polynomial, 
denoted by (p r , is the only irreducible polynomial that divides X r — 1 and does not di¬ 
vide any X s — 1 for s < r. Furthermore, it is known that deg(/ r ) = (p[r) where cp is Eu¬ 
ler's totient function. That is, the cyclotomic degree-d factors of h divide some X r — 1 
with cp(r) < d. In our implementation, we compute the set R = {r < 2d 2 : cp{r ) < d} 
in a naive way, and for each r G R, we exactly compute the r-th cyclotomic polynomial 
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Algorithm 1 Bounded-degree partial factorization of / 

Input: / = E/=i CjX a i with «i < • • • < oc^, and d > 0 ; 
Output: The sets G and H as defined in Eq. (JT]). 


1 
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12 

13 

14 

15 


procedure Partial_factorization(/, d) 

S <— [ciX ai ,.. .,c k X ak ] 

b false 
while |S| > 1 do 
if b then 

8 S cd ( s ) 
if deg(g) > 0 then 

(G,H) PARTIAL_FACTORIZATION(//g, d) 
return (G U {g},H) 

t <— index that minimizes 5 = val(S[f + 1 ]) — deg(S[f]) 

b (d> > 'y(S[f] + S[t + l],d)) 

S[t] i — S[t] + S[t + 1 ] 

S^S\{S[f + l]} 

if b then return ({S[O]}, 0 ) 
else return ( 0 , {S [0]}) 


> Only monomials 


> true or false 


as <p r = (X r - l)/gcd(X r - 1 , riseR,s<r tp s ). Thus, instead of computing the factors 
of gcd(/;, X r — 1 ) as in Lenstra's original algorithm, we compute cp r and simply test 
whether it divides h, using again the fact that <p r divides h if and only if it divides 

mod r 

2.3 Multivariate case 

The multivariate case has the same high-level structure as the univariate case. Namely, 
the computation of the unidimensional factors, that is multivariate factors that can be 
written as f{X) = X a f n {X s ) where X a and are multivariate monomials and f n is 
a univariate polynomial of valuation 0, reduces to the factorization of some univariate 
polynomials, while the computation of the multidimensional factors is based on a Gap 
Theorem. 

The algorithms we implement for these two steps are very close to the one de¬ 
scribed in our previous works |g, 5j. For the multidimensional factors, we apply the 
same modification as for Lenstra's algorithm in the univariate case, that is we first 
write / as a sum of monomials and then iteratively merge the summands. This way, 
we actually build a so-called single-linkage clustering of the set of exponents. 
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3 Examples 


We include an example of timings that is pretty representative of the computation 
times we obtain. It is not comparable to other softwares since in the range of degrees 
that we consider, other software are unable to give any answer. (Often, they even do 
not accept large enough exponents for polynomials.) Also, the nature of the algo¬ 
rithms we use make the computation times very versatile: A very small change in the 
input polynomial may force the algorithm to factor a much larger polynomial, or to 
compute a gcd with much larger inputs. This raises two important questions: Can we 
(automatically) provide an estimation of the size of the polynomials that will appear in 
the computations? Can we modify the algorithms to reduce their versatility, or imple¬ 
ment strategies to compute factors of degree as large as possible while guaranteeing 
reasonable computation times? 

We consider a random polynomial of degree > 1 000 000 and > 10 000 nonzero 
terms constructed as a product of three parts: a product of 5 degree -10 polynomials, 
a product of 3 polynomials of the form X r — 1 with r of order 100000 and a sparse 
polynomial of degree 1000 000 with 40 monomials each multiplied by polynomials of 
degree 20 . Next table gives the timings to compute degree- 5 , 10 , 15 , and 20 factors. We 
observe that the main part of the computation rapidly becomes the gcd computations. 

The code for this example is distributed as part of the Lacunaryx package, in 
the file lacunaryx/bench/bounded_degree_univariate_bench. cpp. The timings were 
obtained on an Intel© Core™ CPU @ 2.60GHz platform with 7.7GB RAM. 


Degree of the factors 

5 

10 

15 

20 

Total time (ms) 

1994 

2924 

12190 

26165 

Non-cyclotomic (ms) 

1948 

2856 

12140 

26061 

Cyclotomic (ms) 

46 

68 

50 

104 

Gcd computations (ms) 

9 i 

749 

9223 

23151 
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